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Abstract. We study the dynamics of homogeneous isotropic FRW cosmologies 
with positive spatial curvature in /(_R)-gravity, paying special attention to the 
existence of Einstein static models and only study forms of f{R) = R" for which 
these static models have been shown to exist. We construct a compact state space 
and identify past and future attractors of the system and recover a previously 
discovered future attractor corresponding to an expanding accelerating model. 
We also discuss the existence of universes which have both a past and future 
bounce, a phenomenon which is absent in General Relativity. 



PACS numbers: 98.80.JK, 04.50.+h, 05.45.-a 

1. Introduction 

After more than one hundred years, Einstein's theory of General Relativity still 
remains our best description of gravity and has up to now survived the scrutiny of 
a multitude of tests, most of which have been on Solar System scales. However, in 
the last decade, high precision astrophysical and cosmological observations appear 
to suggest that General Relativity might be incomplete. In particular, cosmological 
data indicates an underlying cosmic acceleration of the Universe cannot be recast in 
the framework of General Relativity without resorting to additional exotic matter 
components (known as Dark Energy), which have not yet been directly observed. 
Several models have been proposed [1] in order to address this problem and at the 
moment, the one which appears to fit all available observations (Supernovae la [2], 
Cosmic Microwave Background Anisotropics [3] , Large Scale Structure formation [4] , 
baryon oscillations [5], weak lensing [6]), turns out to be the so called Concordance 
Model in which a tiny cosmological constant is present [7] and ordinary matter is 
dominated by Cold Dark Matter (CDM). However, despite its success, the A-CDM 
model is affected by significant fine-tuning problems related to the vacuum energy 
scale, so it is important to investigate other viable theoretical schemes. 

Currently, one of the most popular alternatives to the Concordance Model is based 
on modifications of standard Einstein gravity. Although modifications of Einstein's 
theory of gravity were already proposed in the early years after the publication 
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of General Relativity, a detailed investigation of cosmological models within this 
framework only got underway a few years ago. Such models became popular in the 
1980's because it was shown that they naturally admit a phase of accelerated expansion 
which could be associated with an early universe inflationary phase [8] . The fact that 
the phenomenology of Dark Energy requires the presence of a similar phase (although 
only a late time one) has recently revived interest in these models. In particular, the 
idea that Dark Energy may have a geometrical origin, i.e., a that there is a connection 
between Dark Energy and a non-standard behavior of gravitation on cosmological 
scales is currently a very hot topic of research. 

Among these models the so-called Extended Theory of Gravitation (ETG) and, 
in particular, higher-order theories of gravity (HOG) [9] have provided a number of 
extremely interesting results on both cosmological [10-13] and astrophysical [12, 14] 
scales. These models are based on gravitational actions which are non-linear in 
the Ricci curvature R and/or contain terms involving combinations of derivatives of 
R [15-17]. One of the nice features of these theories is that the field equations can be 
recast in a way that the higher order corrections are written as an energy - momentum 
tensor of geometrical origin describing an "effective" source term on the right hand 
side of the standard Einstein field equations [10,11]. In this Curvature Quintessence 
scenario, the cosmic acceleration can be shown to result from such a new geometrical 
contribution to the cosmic energy density budget, due to higher order corrections of 
the Hilbert-Einstein Lagrangian. 

Unfortunately the analysis of HOG is complicated by the fact that the resulting 
fourth-order field equations are highly non-linear and finding cosmological solutions 
by solving the field equations directly has proved to be extremely difficult. This 
problem has been eased somewhat by using the theory of dynamical systems which 
has over the last few years proved to be a powerful scheme for investigating the physical 
behavior of such theories (see for example [18-25]). In fact, studying cosmologies using 
the dynamical systems approach has the advantage of providing a relatively simple 
method for obtaining exact solutions (even if these only represent the asymptotic 
behavior) and obtain a (qualitative) description of the global dynamics of these models. 
Consequently, such an analysis allows for an efficient preliminary investigation of these 
theories, suggesting which background models deserve further investigation and this 
provides a starting point for the analysis of the growth of structure in HOG [26,27]. 

A remarkable feature of the background dynamics of /(i?) gravity is that there 
are a number of classes of these theories that admit a Friedmann transient matter- 
dominated decelerated expansion phase, followed by one with an accelerated expansion 
rate [18,28]. The first phase provides a setting during which structure formation can 
take place and this is followed by a smooth transition to a Dark Energy like era which 
drives the cosmological acceleration. It would therefore be of great interest if orbits 
could be found in the phase-space of f{R) models that connect such a Friedmann 
matter dominated phase to an accelerating phase via an Einstein Static solution in 
a way which is indistinguishable (at least at the level of the background dynamics) 
to what occurs in the LamfodaCDM model of General Relativity [29]. If a similar 
evolution exists in an /(i?) model, we would expect the associated Einstein Static 
solution to be a saddle point (as it is in General Relativity) and consequently unstable. 
It therefore follows that a careful examination of the existence and stability of the 
Einstein Static model in the more general setting of modified gravity is of critical 
importance in determining whether ACDM - like cosmologies are a generic feature of 
such gravitational theories. 
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In General Relativity, the issue of stability of the Einstein Static model has been 
studied several times since the classic paper by Eddington [30] in the 1930s, where it 
was shown that such models are unstable with respect to homogeneous and isotropic 
perturbations, exactly the feature that allows a transition between a decelerated 
expansion era to one which is accelerating. However, later work by Harrison [31] and 
Gibbons [32], which extended Eddington's work, considered generic inhomogeneous 
and anisotropic perturbations of an Einstein Static model filled with a perfect fluid 
and found that provided the sound speed satisfies Cg > , these models are marginally 
stable with respect to such perturbations. The reason for this "non-Newtonian" 
stability stems from the fact that in General Relativity, the Einstein Static universe 
is spatially closed, and therefore has a maximum scale associated with it, which is 
greater than the largest physical scale of the perturbations. Since the Jean's length is 
a significant fraction of this maximum scale, perturbations in the fluid oscillate, rather 
than grow, leading to the conclusion that in General Relativity at least, the Einstein 
Static solution is marginally stable with respect to such perturbations. This result 
was also found more recently by Barrow et. al. [33] 

Recently [35] the existence and stability of the Einstein Static models in the more 
general setting of f{R) gravity was examined. It was found that only one class of 
f{R) theories admits an Einstein Static model, and that this class is neutrally stable 
with respect to vector and tensor perturbations for all equations of state on all scales. 
Scalar perturbations are only stable on all scales if the matter fluid equation of state 
satisfies Cg > ~ 0.21. This result is remarkably similar to the GR case discussed 

above. 

In this paper we use the theory of dynamical systems to examine in detail the 
background evolution of this special class of f{R) theories. This analysis compliments 
the work in [35] and provides a more complete picture of how the Einstein Static 
solution fits into the overall structure of the phase space of solutions in f{R) gravity. 

The following conventions will be used in this paper: the metric signature is 
( — h H — h); Latin indices run from to 3; units are used in which c = SttG = 1. 



2. The Field Equations 



For homogeneous and isotropic spacetimes, a general form of the action for fourth- 
order gravity is given by 



A 



-g [fiR) - 2A + 



(1) 



where £,„ is the Lagrangian of the matter fields. The fourth-order field equations can 
be obtained by varying (jT]): 

A T™ 

Gab + 9abjj = -p- + T^b = T^b > (2) 

where primes denote derivatives with respect to R. Here is the total effective 
energy momentum tensor composed of the ordinary matter energy momentum tensor 
and the correction term (often referred to as the "curvature fluid"): 



= — 

ab J, 



i^gabU^Rn + fUalai-g^'gab) 



(3) 
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As shown in [35], if we wish to keep the extra degrees of freedom in fourth-order 
gravity, the existence of an Einstein static universe imposes a strong constraint on the 
function f{R), which was found to be 

/(i?) = 2A + /Ci?^(i+"') , w 7^ -1. (4) 

Here we assume that the standard matter is a barotropic perfect fluid with equation 
of state 

P"' = wp"' , (5) 

where and are the standard matter density and pressure respectively. The 
effective field equations will look like the ones for i?"-gravity, but subject to the 
constraint 

n^'^il + w), (6) 

which guarantees the existence of an Einstein Static solution. This means the equation 
of state is fixed as a function of n, or in other words: for a given value of n we can 
only find Einstein static solutions if the equation of state satisfies (|6]). As we know, in 
order to satisfy all the energy conditions for standard matter, w has to take values in 
[—1/3,1]. We will therefore only consider the range n S [1,3], and in particular will 
analyze in detail the cases of dust (n = 3/2, w = 0) and radiation (n = 2, w = 1/3). 
In what follows we will use equation ^ to eliminate w from the field equations. 

For the homogeneous and isotropic spacetimes, the independent field equations 
for i?"-gravity are (with n given by (6): 

• The Raychaudhuri equation 



e + ie^--i?-(.-i)-e + ^ = o, (7) 



where Q is the volume expansion which defines a scale factor a(t) along the fluid 
flow lines via the standard relation O = 3d/ a. 

• The Friedmann equation 

where R = (6K)/a^ is the curvature of the 3-spaces and k G (1,0,-1) denotes 
closed, flat and open universes respectively. 

• The trace equation 

i?_lJn-2^ . R^ ,2 (2-n) 

• The conservation equation for standard matter and the propagation equation for 
the 3-curvature are 

p'" = -|v"e , R^-^m (10) 

respectively. 

Combining the Raychaudhuri and Friedman equations, we get 

i? = 29 + ^62+ J?, (11) 
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which is equivalent to the general definition of the Ricci scalar in terms of the scale 
factor in a FRW spacetimes, given by 

„ f ci a? k\ , ^ 

i? = 6 - + - + - . (12) 

For our analysis of the dynamical system, which will be presented in the next section, 
it is useful to complete the square of equation (8) and rewrite the Friedmann equation 
as 

y 2 Rj 2 4 R^ 2n nR''-^ ^ ' 

We can easily see that all the additive terms on both sides of the equation above 
are strictly non-negative for all spacetimes with non-negative Ricci scalar R and 3- 
curvature 7?. This will help us to compactify the full FRW state space for non-negative 
spatial curvature and Ricci scalar, and thus allow us to analyze the dynamics of the 
system in a compact framework. 



3. Dynamics of non-negative curvature FRW spacetime 

We now study the dynamics of FRW models with non-negative spatial curvature 
and Ricci scalar, which include the Einstein static universe. In order to convert the 
equations above into a system of autonomous first order differential equations, we 
define the following set of normalized variables 0; 

-(n-1), y= ^(71-1), z= — ^, (14) 



2RD^ " 2nD^^ nR'^'^D^ 

3R n ^ 



together with the normalized time variable 
, _ d _ 1 d 
^ dr ^ 15 dt ■ 
Here the normalization D is chosen as 



(15) 



D = 



\ 



^ 3(n-l)i?\^ 3~ 



in order to compactify the variables as shown below. In terms of these variables, the 
Friedmann equation (|13p becomes 

x^+y + z = l, y,z>0, (17) 

and from the definition of the normalization we get another constraint 

{Q + xf+K=l , Jv >0. (18) 

The above constraints show that the variables must be compact and take the following 
ranges: 

X e [-1, 1] , 2/ e [0, 1] z e [0, 1] , Q e [-2, 2] , if e [0, 1] . (19) 

f It is important to note that tliis choice of variables will exclude General Relativity, i.e., the case of 
n = 1. See [29] for the dynamical systems analysis of the corresponding cosmologies in GR. 
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Using the five compact variables together with the two constraints, we reduce the 
complete dynamical system to a 3-dimensional one, and the cosmological equations 
become equivalent to the autonomous system 

Q' = [(3 - n)x^ _ n(y _ 1) _ 1] ^ + [(3 - - n{y - 1) 



1 

3 ^ 



-1 



ny 



1 



g (3-n)(x + 0) + ^ 



(n^ - 2n + 2) 



n — 1 



ny 



1 



Qx 



';Qny{l - y) 



(20) 



= ^(3 'n){Q + x) + '^ [n(2 ~ y) - b] 



3 

Qx 



^ - 2/) - 3] + - 



n(n- 2) 



- ?i + 2 



4. Qualitative Analysis of the state space 

The complete description of all the equilibrium points of the dynamical system 
described by the above equations is given in Table I. We note that the points 
A and B correspond to the similarly labeled points in [36]. The line CC has 
been labeled as such since it includes the point C found in [36] in the flat limit 
+ = \{2 — n)Q\ — > 1. Strictly speaking the system (|^D|) has an additional 
equilibrium point, but its coordinates only lie in the state space for n > 3, which for 
us corresponds to an unrealistic equation of state w > 1. Note that the subscript e 
only labels expanding (e = +1) and collapsing (e = —1) models for the point A. Since 
the other points correspond to asymptotically Minkowski spaces, the subscript merely 
labels the different limits from which these asymptotic solutions were obtained. The 
line CC has an expanding branch CC^^^ with Q > 0, a collapsing branch /^C^"" with 
Q < 0, and in between the Einstein Static model £5, with Q ~{). The stability of the 
equilibrium points and the line is summarized in Table II. We will now briefly discuss 
several interesting features of these equilibrium points. 

(i) The position of the equilibrium points Me in the state space does not depend on 
the equation of state for the standard matter. In fact, the existence of these points 
is a generic feature of i?" gravity, and they form a global source (e = +1) or sink 
(e = — 1) in the state space considered here. This implies that for any equation of 
state, there always exists the possibility of the Universe originating from a vacuum 
Minkowski model in the infinite past, or evolving towards a vacuum Minkowski 
model in the infinite future. 

(ii) The position of the equilibrium points also does not depend on the equation 
of state for standard matter. However, the stability properties of these points 
do depend on the equation of state. For matter with high negative pressure 
(—1/3 < w < —1/6), or for very stiff matter (2/3 < w < 1), these points are 
unstable saddles. For other values of the barotropic index, these points represent 
a global source/sink depending on the sign of e. 

(iii) The points only exist for the specific range of the barotropic index w G 
[—1/3, 2/3], and they are always unstable. 
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(iv) The points A^^ exist for w G [—1/6, 1] only, and they have the interesting property 
of exhibiting de Sitter (dS)/anti-de Sitter (AdS) hke acceleration/deceleration 
depending on the barotropic index. In this case however, the acceleration has a 
power law behavior instead of the exponential behavior of dS/AdS. These points 
are vacuum solutions, and the expanding point is a global sink if n > P+ w 1.37 
01 w > w+ w —0.09. These points are also a generic feature of i?" gravity, and 
the existence of these points implies that even with an ultra-violet correction to 
General Relativity, one can evolve towards late time acceleration for the ranges 
of w given in Table I. 

(v) The line CC is entirely an artifact of the n — w correspondence (j6|) , and is absent 
in a general i?"-gravity. All the points on this line are expanding/collapsing non- 
accelerating solutions, containing the Einstein static solution as the only static 
point. For n G [1, P+] ~ [1, 1.37] (see table caption), the expanding branch of the 
line forms spiral sinks, while the collapsing branch of the line forms spiral sources. 
The Einstein static point in between these branches forms a center, that means 
it is neutrally stable for this range of n. For n £ [P+,3/2] a bifurcation appears 
on the line: for Q G [—Qb^Qb], the expanding (contracting) points on the line 
remain spiral sinks (sources), and the Einstein static point remains a center. The 
points with |Q| > Qb however now are saddles. Here Qb is the bifurcation value 
defined by 



We can see that the bifurcation enters the state space from both endpoints of the 
lines for n = P_|_. As n increases, Qi, and —Qb move closer towards the center of 
the line, and for n = 3/2 the two bifurcations iQf, merge at Q = 0. For equations 
of state stiffer than dust (n > 3/2), all the points on the line are saddles, including 
the Einstein static model. 

(vi) The state space in general contains several past and future attractors. 
Interestingly, we find that for any equation of state, there is no expanding past 
attractor. The only possible past attractors are the asymptotically Minkowski 
points and Af+, and the decelerating collapsing point For w < 0, the 

points on the collapsing non-decelerating line £C'^°" may also be past attractors. 
This means that in this class of models, there is no Big Bang scenario for any 
realistic equation of state w G [—1/3,1], but only possible bounce scenarios or 
expansion after an initial asymptotic Minkowski phase. Figure I shows the state 
space for dust {n = 3/2), with all the equilibrium points. There are bounce 
scenarios along those trajectories where the quantity Q changes sign. Figure 4 
shows a solution which is bouncing in past as well as in the future. These kind 
of solutions are not present in General Relativity. 

Stability of the Einstein Static model 

We will now discuss the stability of the Einstein static point in more detail. The 
stability of the Einstein Static model against general covariant, gauge-invariant linear 
perturbations was studied in [35] . The standard procedure of harmonic decomposition 
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Figure 1. The complete state space for dust with all the equilibrium points. We 
can see the trajectories coming out from the global sources C+, Af+ and A.- and 
going to the sinks Af- and . Along the trajectory from A.— to the 
quantity Q changes sign twice which shows bounce scenarios in past and future. 



is used, employing the trace-free symmetric tensor eigenfunctions Q of the spatial 
Laplace-Beltrami operator: 

V^O = -^Q , Q = , (22) 
% 

where the constant oq defines the length scale for Einstein Static universe. For the 
spatially closed models, the spectrum of eigenvalues is discrete and given by [31,37] 

fc2 = n(n + 2), (23) 

where the co-moving wavenumber takes the values h — 1,2,3, We note that 
the dynamical systems analysis discussed in the present paper corresponds to the 
homogeneous mode with comoving wavenumber ri = in the harmonically decomposed 
perturbation equations. This mode was not studied in [35], since it corresponds to a 
change in the background of the linear perturbation theory. 

It can be easily seen from equation (53) in [35] that the scalar perturbations for 
n = have a growing and a decaying mode if w > 0. If w < on the other hand, 
the density gradient has two purely oscillatory modes. These results agree with the 
stability of the Einstein Static model in the dynamical systems context analyzed here, 
where we found that the Einstein Static model is a saddle for w > 0. For w < 0, since 
all the points of the line CC at one side of Einstein Static are sources and other side of 
Einstein Static are sinks, the point Einstein Static becomes quasiperiodic in nature, 
in the sense that orbits emerge from a source, spiral around Einstein Static for some 
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x(1) 




Figure 2. The complete state space for matter with w = —1/9, with all the 
equilibrium points. The expanding part of line CC now behaves as a spiral sink. 
Also the points A.— and arc now saddles. 



time and goes to a sink. Figure 3 shows such an orbit, where the initial conditions 
are taken very close to the Einstein Static point on the plane perpendicular to the 
line CC containing the Einstein Static point. This quasiperiodic nature reflects in the 
oscillatory modes of the scalar perturbations in [35]. 

As we can easily see from the dynamical system equations, for dust-like matter 
{w = 0) the £S point is a double bifurcation point, and consequently the linear theory 
breaks down. The physical interpretation is as follows: from equation (53) in [35] we 
can easily see that for w = and n = 0, the linear perturbation is independent of 
time, i.e., it is a constant. Hence to study the actual behavior of the homogeneous 
perturbation for dust, we must consider higher order corrections. 

4-2. The state space for radiation 

For n = 2 (or w = 1/3) the state space has some interesting properties. Though 
qualitatively the nature of all equilibrium points arc similar to that of dust, due to 
vanishing of the trace of the energy momentum tensor for standard matter, the surface 
defined by a; = becomes an invariant subspace, containing the global sink/source 
Ae, global saddles and the Einstein static point £S which is also a saddle. By 
definition, x = corresponds to i? = 0, and therefore the points Ae are not accelerated 
expansion/collapse solutions but they are instead static Minkowski points. Figure 5 
shows the phase portrait of this subspace. It is interesting to note that in this subspace 
the points Be behaves like source/sink, however in the complete state space they are 
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x(1) 




Q(t) 



Figure 3. The quasiperiodic nature of the Einstein Static point for w = —1/9. 
Here the initial conditions arc taken very close to the Einstein Static point on the 
plane perpendicular to the Unc CC containing the Einstein Static point. 



It 1n 




Oil) Q(t) 

Figure 4. Example of bouncing solutions for dust (left) and radiation (right). 
The dust solution shows a bounce in both past and future which is not possible 
in General Relativity. 



saddles. All the trajectories from A- to A+, and from B+ to B- admit bouncing 
universes. Also figure 4 shows a solution for radiation like matter linking to £_ 
which has a bounce in the future. 
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Table 1. Equilibrium points of the FRW state space. We have abbreviated 
'^W = 3(n-?K2n-i) ■ = ('^ " 3)('^ " 1)' '^"d P+ = (1 + V3)/2 1.37. The 
quantities ao, 01,02,13 are constants of integration. 



Point {Q,x,y) 


constraints 


Solution/Description 


K (0, e, 0) 


n e [1,3] 


Vacuum Minkowski 


J^e (2e, -e, 0) 


r~ ri Ql 

n t [i, oj 


Vacuum Minkowski 


B, ( ^ e^, 0] 

^ yS—n n— 3' J 


n G [1,2.5] 


Vacuum Minkowski 


A ( ^ 2n-l , n-2 8n^-14ri+5^ 
\^''3(n-l) ' ^3(n-l) ' 9(n-l)^ ^ 




n e [1.25,3] 


Vacuum, Flat, AccelerationT^ 
Decelerating for P+ < n < 2 
a{t)=ao {ai + k{n)ty^''^''^ 


Line 

CC (0, Q(n 1), 




IQI < 2^ forne [1,P+] 

v^^Ufo'^-^t^+'S] 


Non- Accelerating curved 
ait) = flat + 03 , > 



Table 2. This table summarizes the nature of the equilibrium points. As in the 
previous table, we have abbreviated P+ = (1 + V3)/2 1.37. The physically 
interesting values correspond to n = 3/2 (dust) and n = 2 (radiation). Note for 
n = 3/2, the Einstein Static point becomes a bifurcation, and its stability cannot 
be determined within the linear theory (see text). 



point 


type 


range of n 


(1,5/4) (5/4, P+) (P+,3/2) (3/2,5/2) (5/2,3) 


A+ 


expanding 




saddle 


sink 


A- 


collapsing 




saddle 


source 


B± 


static 


saddle 




c+ 


static 


saddle 


source 


saddle 


C- 


static 


saddle 


sink 


saddle 


M+ 


static 


source 


N- 


static 


sink 




expanding 


sink 


sink (for Q < Qt) 
saddle (for Q > Qi,) 


saddle 




collapsing 


source 


source (for \Q\ < Qi,) 
saddle (for |Q| > Qt,) 


saddle 
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Figure 5. The x = invariant subspace for radiation, containing the points 

Be and £S. It is interesting to note that in this subspace the points Be behaves 

hke source/sink, however in the complete state space they are saddles. 



5. Discussion and Conclusions 

In this work, we have studied the state space of the class of isotropic FRW models in 
/(i?)-gravity. We were particularly interested in the stability of the Einstein Static 
model, which in general only exists for the specific form of f{R) = i?", with n 
constrained to be a function of the equation of state parameter w as shown in [35]. 
For this reason, we here considered i?"-gravity with the constraint n = 3/2(1 + w) 
only. 

We have identified the past and future attractors of the state space of these 
models. We have found that the Einstein Static model is an unstable saddle point 
for all equations of state stiffer than dust (w > or n > 3/2), and a neutrally stable 
center for w G (—1/3,0), i.e. n £ (1,3/2). This is different to General Relativity, 
where the Einstein Static model is an unstable saddle for all w g (—1/3, 1). We have 
also numerically found bouncing orbits that link the decelerating collapsing point A- 
to the accelerating expanding point A-^- via an asymptotically Einstein static phase 
represented by the saddle point £S. 

We note that the stability of the Einstein Static model obtained here is in agree- 
ment with the analysis of inhomogcncous perturbations [35], provided we carefully 
examine the relationship between the dynamical systems analysis here and the har- 
monically decomposed linear perturbations described in [35] . The homogeneous mode 
with comoving wavcnumbcr h = was excluded in [35], since it corresponds to a 
change in the background of the linear perturbation theory. It is this mode however 
that corresponds to the stability properties of the dynamical system examined here. 
This mode is neutrally stable if w < 0, and unstable with a decaying and a growing 
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mode ii w > 0. These results exactly match the results from the dynamical systems 
analysis presented in this paper. 
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